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ABSTRACT 

It is shown here that in a flat, cold dark matter (COM) dominated Universe with positive 
cosmological constant (A), modelled in terms of a Newtonian and collisionless fluid, particle 
trajectories are analytical in time (representable by a convergent Taylor series) until at least a 
finite time after decoupling. The time variable used for this statement is the cosmic scale fac¬ 
tor, i.e., the “a-time”, and not the cosmic time. For this, a Lagrangian-coordinates formulation 
of the Euler-Poisson equations is employed, originally used by Cauchy for 3-D incompress¬ 
ible flow. Temporal analyticity for ACDM is found to be a consequence of novel explicit 
all-order recursion relations for the a-time Taylor coefficients of the Lagrangian displacement 
field, from which we derive the convergence of the a-time Taylor series. A lower bound for 
the a-time where analyticity is guaranteed and shell-crossing is ruled out is obtained, whose 
value depends only on A and on the initial spatial smoothness of the density field. The largest 
time interval is achieved when A vanishes, i.e., for an Einstein-de Sitter universe. Analyticity 
holds also if, instead of the a-time, one uses the linear structure growth D-time, but no simple 
recursion relations are then obtained. The analyticity result also holds when a curvature term 
is included in the Eriedmann equation for the background, but inclusion of a radiation term 
arising from the primordial era spoils analyticity. 
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1 INTRODUCTION 

Recent observational results (Planck Collaboration XVI 2014) indicate that we live in a spatially flat ACDM Universe which is nowadays 
dominated by a cold dark matter (CDM) component and a cosmological constant, A > 0. In a primordial era, matter was tightly coupled to 
radiation via electroweak interactions. This tight coupling prevented matter to cluster significantly at early times since the radiation pressure 
acted as a counter force to the gravitational force. As the Universe expands its mean energy density decreases, which eventually lead to 
a freeze-out of these electroweak interactions. This freeze-out of the interactions (non-instantaneous in space and time) happened roughly 
380,000 years after the Big Bang, eventually enabling the matter to cluster. This epoch, usually called decoupling or recombination, marks 
the beginning of cosmological structure formation with initially smooth matter density fluctuations. 

The observed large-scale structure of the Universe is mainly the result of the gravitational instability. To study the structure formation, 
it is usually assumed that the CDM component behaves as a Newtonian pressureless and curl-free fluid (so-called cosmological dust). Such a 
fluid is governed by the Euler-Poisson equations, together with the Friedmann equation, where the latter describes the background evolution 
of the expanding Universe. Since that background evolution is by definition exactly homogeneous and isotropic, it is parametrized by a single 
function (of cosmic time t): the cosmic scale factor a{t). 

Analytical models for cosmological structure formation can be formulated in either the Eulerian or Lagrangian frame of reference (for 
some reviews, see e.g. Bemardeau et al. 2002; Bemardeau 2013). The class of analytical models of the latter is dubbed the Lagrangian 
perturbation theory (LPT) (Buchert & Goetz 1987; Buchert 1989; Buchert 1992; Bouchet et al. 1992; Bildhauer et al. 1992, Bouchet et 
al. 1995; Catelan 1995; Ehlers & Buchert 1997). When using this technique, the only dynamical variable is the displacement, which represents 
the gravitationally induced deviation of the particle trajectory field from the homogeneous background evolution of the Universe. To solve 
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the Euler-Poisson equations in Lagrangian space, one usually requires some power series Ansatz for the displacement (for more information 
see also below in the Introduction). 

The LPT is widely applied in cosmology. Buchert (1992) showed that the ZeTdovich approximation (ZeTdovich 1970) is actually an 
instance of the first-order LPT (see also Novikov 1970). Simply and elegantly, it quite well describes the gravitational evolution of CDM inho¬ 
mogeneities, as long as the trajectories of fluid particles do not cross. A Lagrangian approach has been also used in the so-called cosmological 
reconstruction problem in an expanding Universe (Frisch etal. 2002; Brenier etal. 2003), where it is shown that, in so far as the dynamics are 
governed by the Euler-Poisson equations, the knowledge of both the present highly non-uniform distribution of mass and of its primordial 
quasi-uniform distribution, uniquely determines the inverse Lagrangian map, defined as the transformation from present Eulerian positions 
to their respective initial positions. The LPT is also an important tool in numerical investigations. In recent years, second-order LPT (2LPT) 
has been successfully used to set up initial conditions for numerical A-body simulations (see e.g. Crocce, Pueblas & Scoccimarro 2006). 
The so-called biasing model, whose task is to find a relationship between the visible matter and the dark matter distribution, yield good 
predictions to cosmological observations when formulated in Lagrangian space (Mo & White 1996; Matsubara 2008b). Also, semi-analytical 
approaches to LPT deliver statistical estimators, such as the matter power spectrum and bispectrum, which compare favourably with results 
from A-body simulations (e.g., Matsubara 2008a; Rampf & Wong 2012; Vlah, Seljak & Baldauf 2014). These approaches make use of the 
LPT up to the fourth order (Rampf & Buchert 2012; Rampf 2012). Most of the LPT calculations are performed for a flat CDM universe with 
vanishing cosmological constant, i.e., for an Einstein-de Sitter (EdS) universe, and only few investigations, up to the third order, are known 
for a ACDM Universe (see, e.g., Matsubara 1995; Lee 2014). The LPT has also been generalised to General Relativity, for a non-perturbative 
formulation see Buchert & Ostermann (2012) and Buchert, Nayet & Wiegand (2013). For a general relativistic treatment within the ACDM 
model see Rampf & Rigopoulos (2013), Rampf (2013), Rampf & Wiegand (2014). 

Although LPT is widely applied in cosmology, very little is known about its convergence properties, which requires some knowledge 
about terms of arbitrarily high order. There are of course some exceptions, for example, Sahni & Shandarin (1996) investigated the case of 
an initial “top-hat underdensity” that is initially discontinuous, and found that low-order perturbation worked better than a higher-order one, 
which they regarded as a possible evidence for “semiconvergence” of the perturbation series (see also Moutarde et al. 1991). Recently, a 
novel Lagrangian-coordinates approach has been used to show that particle trajectories for an EdS universe are time analytical until a finite 
time (Zheligovsky & Frisch 2014). Their approach is based on a little-known Lagrangian formulation of ideal fluid flow derived by Cauchy 
in 1815 (Frisch & Villone 2014). In this paper we show that the approach of Zheligovsky & Frisch (2014) can be extended to a ACDM 
Universe (and even beyond, see below), derive novel recursion relations and prove time analyticity in the cosmic scale factor time, here used 
both as a time variable and as an expansion parameter. 

In the LPT literature, analytical results for ACDM usually employ an expansion in powers of small displacements, which amounts to 
performing a (Taylor) expansion in powers of the initial peculiar gravitational potential, ~ 10“®. By this technique one seeks (rea¬ 

sonably well) approximate solutions of the non-linear Euler-Poisson equations. As a consequence of such an approximation technique, one 
is forced to solve at each order in perturbation theory a second-order partial differential equation for the time coefficient of the displacement 
(the fastest growing solutions of these time coefficients are usually denoted with D(t), E{t), etc.). Because of that it is impossible to obtain 
recursion relations for ACDM by the use of such a “small displacement” expansion. In this paper, by contrast, we perform a time-Taylor 
expansion and seek for exact analytic solutions for the fully non-linear Euler-Poisson equations in Lagrangian space. As a consequence of 
our expansion scheme, the displacement is represented in terms of a power series in the cosmic scale factor (even for a ACDM Universe!), 
and there is thus no need to explicitly solve partial differential equations for higher-order time coefficients. 

This paper is organised as follows. In Section|2]we present various forms of the Euler-Poisson equations in a Eulerian-coordinate system 
and we show that regularity of the solution at short times requires certain slaving constraints on the initial conditions. In Section[3]we turn 
to the Euler-Poisson equations in Lagrangian coordinates, which take a particularly simple form when using the cosmic scale factor as the 
time parameter. Recursion relations are derived in Section |4l from which time-analyticity is derived in Section Further results related to 
time-analyticity are discussed in Section|^ This includes the dependence of the time of guaranteed analyticity on the cosmological constant 
A (Section [6. lb . the absence of shell-crossing during the time interval where analyticity can be proved (Section [6.2b . the analyticity in the 
linear growth time variable D (Section (63) and the persistence of analyticity when curvature effects are included, as well as the problem 
that arises with radiation effects (Section r6.4b . Finally, in SectionlT] we summarise our main findings and highlight various challenging open 
problems, such as using analyticity to develop a semi-Lagrangian numerical approach to the cosmological reconstruction. 


2 THE VARIOUS FORMS OF THE EULER-POISSON EQUATIONS AND THE SLAVING CONDITIONS 

A flat, matter dominated Universe may be studied as a Newtonian and collisionless fluid whose governing equations are the Euler-Poisson 
equations. For a flat Universe with cosmological constant, the Euler-Poisson equations are (Peebles 1980) 


dtU+(tj-Vr.)u = -Vr^g, 

(la) 

dtg + Vt- • = 0, 

(lb) 

= A-kGq — 3A . 

(Ic) 
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Here, r is the proper space coordinate, t the cosmic time, U the velocity of the fluid, its density q\ the gravitational potential is denoted by 
(j}g. As to the cosmological constant, it is here denoted by 3A. Observe that, as long as we use the cosmic time as independent time variable, 
the various dependent variables are surmounted by a tilde, which will be dropped when we change the time variable from cosmic time to 
cosmic scale factor a. As usual, we decompose the mass density, fluid velocity and gravitational field respectively in a sum of two terms, one 
describing the effect of the uniform background expansion, the other the fluctuations against its background. 


Q = qU) [l + ^] > 


U = —r + au , 
a 


— 4>g + tpg ■ 


( 2 ) 


Here, 5 is the density contrast and aii the proper peculiar velocity. Furthermore, a{t) is the cosmic scale factor which parametrizes the global 
background expansion governed by the (first) Friedmann equation 
2 

(3) 


— a ^ + A, 


where we have used the fact that the background mass density is git) ~ a ®(f), and we have set StvGqo /3 = 1 for convenience. The solution 
for the cosmic scale factor is easily obtained from the Friedmann equation, namely 

a{t) = A“^^®[sinh(3\/Af/2)]^^® , (4) 

which is proportional to for small t (see e.g., Kofman, Gnedin & Bahcall 1993). 

After performing the aforementioned decomposition, derived in Appendix lAl we obtain the Euler-Poisson equations for the (comoving) 
peculiar velocity n 


OtU I'll ■ ^ 7 x) 'll — ‘2 ' H ^ 7 x^g ■ 

a 


dtS + Va 


(1 + 5){ij 


= 0 , 


(5a) 

(5b) 

(5c) 


-72 _ 3 ; 

Note that the peculiar Euler-Poisson equations depend on the comoving coordinate x = r/a. Although these equations are indeed valid for a 
fluid model with cosmological constant, the latter does not explicitly appear in the peculiar approach, but only implicitly through Friedmann’s 
background evolution equation 13. 

Before considering the Euler-Poisson system in Lagrangian space, let us briefly discuss the linearised system in Eulerian space. 

Formally linearising around the steady state n = 0, <5 = 0, we obtain the single differential equation for the density contrast, 

( 6 ) 

a 2a-^ 

The solution of is most easily obtained by changing the time variable from cosmic time t to the cosmic scale factor a, here called a-time. 
We thus set S{t) = 5{a{t)). We then have dt = ada- The above differential equation is then 


(1 + Aa^) 


da 


2a 




(7) 


where we have used the Friedmann equation This equation has two solutions. One is called the decaying mode; it behaves as a for 
small a and thus blows up when a —>■ 0, thereby invalidating the linearisation. The other one, the growing mode, here called the linear growth 
function, is taken to be D{a) — a^/1 + Aa^ 2 T 1 (3/2, 5/6,11/6, —Aa®), where 2 F 1 is the Gauss hypergeometric function (Demianski 
et al. 2005; Enqvist & Rigopoulos 2010; Hamber & Toriumi 2010). This solution is analytic around a = 0, has the small a expansion 
D{a) = a — (2/11) + h.o.t. and thus essentially reduces to the cosmic scale factor at short times. This is the physically appropriate 

solution, in which density fluctuations are growing linearly with a. 

However, when studying the well-posedness of the non-linear Euler-Poisson system dSaHSct and, subsequently, when looking for 
analytic solutions in Lagrangian coordinates, a linear approach ceases to be appropriate. The latter just gives us a signal that the appropriate 
time variable might be proportional to or to a or D at short times. The non-linear theory can be developed in all three such time variables 
but, for our purposes, it takes its simplest form when using the o-time (and not just at short ones). Indeed the Euler-Poisson system ( l53f53 
can be rewritten in the following form: 


(1 + Aa®) [daV + {v ■ Vx) v] = —^ (v -f VxV^s) “ SAa'^v , 

daS + Vx • [(1 + 5) tt] = 0 , 



where we used a new velocity variable v, the a-derivative of comoving particle positions, related to the old velocity by 


( 8 a) 

( 8 b) 

( 8 c) 


'u{x,t) = dv(x,a(t)), (9) 

and we have set ifig{x,t) = {3/2)(pg{x, a{t)). Note that, with respect to A, the situation in fSalf^ is the opposite of what we had 
with ([SaJlScJ: A appears now explicitly, but the Friedmann equation can now be omitted as long as we do not want to revert to the cos¬ 
mic time variable. 
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An important feature of i fSaHS^ is the presence of cosmic scale factors in the denominators of the r.h.s. of and It indicates that 
the solution will be singular at a = 0 unless the corresponding numerators also vanish, that is, we have to satisfy the two following slaving 
conditions (Brenier 1987) at o = 0 (denoted by the superscript (init)); 

^(init) ^ ^ . (10) 

The second slaving condition immediately implies the curl-free character of the initial velocity, which persists by in Eulerian (but not in 
Lagrangian) coordinates. Traditionally, this irrotationality is often derived by linearising the Euler-Poisson equations and showing that the 
rotational component of the velocity decays in time (see e.g., Bernardeau et al. 2002). This linearisation would be justified if the velocity v 
would be small at short times as indeed the density contrast 5 is. However, the velocity v and the gravitational force do not vanish as 

a —>■ 0. Thus, a linearisation-based argument is questionable. 

The cosmic scale factor a, being dimensionless, is conveniently normalised to unity at the present-time epoch. At decoupling it then has 
a value of about 10“^. Here, we let the cosmic scale factor start at the value zero, while using the Euler-Poisson model with slaving. As a 
consequence, the whole primordial (pre-decoupling) cosmology is just reduced to the prescription of the slaving conditions Oni). We shall 
come back to this issue in the Concluding Remarks |7] 


3 THE LAGRANGIAN FORMULATION OF ACDM 


Time-Taylor expansions can be carried out in either Eulerian or in Lagrangian coordinates. For such expansions to be convergent, i.e., for 
the radius of convergence not to vanish, it is much preferable to work in Lagrangian coordinates. Indeed, if the initial conditions have spatial 
derivatives only up to some finite order, the Eulerian solutions will have time derivatives only up to the same order. The reason is that, when 
one takes a time derivative in Eulerian coordinates, this generates a space derivative, because of the [v ■ Vx)i> term. Hence time-Taylor 
coefficients will not exist beyond that order. This does not happen in Lagrangian coordinates, where a limited amount of smoothness in 
the initial data suffices to ensure time-analyticity (Zheligovsky & Erisch 2014). Even when the initial data are spatially analytic, so that 
the solutions are time-analytic in both Eulerian and Lagrangian coordinates, the Lagrangian radius of convergence, being controlled by the 
largest strain in the initial data, is typically much larger than the Eulerian one, which depends on the largest initial velocity present (for more 
details on such matters, see Podvigina et al. 2015). 

We now turn to the Lagrangian formulation of the Euler-Poisson equations fSaH^ . For this we use the (direct) Lagrangian map 
q I—>■ x{q, a) from the initial (a = 0) position q to the Eulerian position x at time a. The velocity v, as it is apparent from its definition 
is simply the Lagrangian a-time derivative of the position; more precisely, v{x{q, a), a) = dax{q, a), where da is the Lagrangian a-time 
derivative. At initial time, a = 0, the velocity is 

= n(a:(< 3 r, 0 ), 0 ) , (11) 


which coincides with the Eulerian one. We shall also use the Jacobian matrix {q, a), where denotes the Lagrangian space derivative 
with respect to qj, and its determinant J = det (V^Xi), the Jacobian. Note that after shell-crossing, i.e., the first vanishing of J, the map 
ceases to have a unique inverse. 

As an intermediate step to a fully Lagrangian formulation of the Euler-Poisson equations, we rewrite the momentum equation 
and the Poisson equation in a mixed form: some terms are more naturally expressed in Lagrangian coordinates and others in Eulerian 
coordinates; composition with the direct or inverse Lagrangian maps is here understood where appropriate; 


(1 + Aa®) daaX = {daX + Va:‘PgJ 


— SAa^dl^x . 


{1 + S)J = 1, 
<5 




(12a) 

(12b) 

(12c) 


Equation dl2bl l is a straightforward expression of mass conservation in Lagrangian coordinates, which need not be derived from its differential 
Eulerian counterpart. It is however invalid after shell-crossing. Due to the appearance of several velocities (multi-streaming), there are several 
branches of the inverse Lagrangian map. In that case dl2b> must be replaced by 


(1 + 5) = E 

n 



(13) 


where n labels the various Lagrangian locations having the same Eulerian location (Buchert 1995). It is likely that each branch of the 
multi-streaming system still satisfies the Euler-Poisson equations with the density given by GS, but this requires further investigation. 

The derivation of the full Lagrangian formulation of the Euler-Poisson equations has two parts. The first part is almost identical to 
Cauchy’s 1815 derivation of the Cauchy invariants equation for incompressible fluids (see Frisch & Villone 2014). Since it is already given 
in Zheligovsky & Frisch (2014) for the case of an EdS universe, we shall just recall the ideas briefly and give the final Cauchy equation. We 
observe that in the r.h.s. of dl2ab . because of the curl-free character of the velocity d^x, all the terms are Eulerian gradients. Multiplication 
by the Jacobian matrix changes these Eulerian gradients to Lagrangian gradients. These are then made to vanish by taking a Lagrangian 
curl. After this, a quadratically non-linear equation is obtained for the Lagrangian map, which involves first and second time derivatives. 
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Cauchy observed that this equation can be exactly integrated in time. This “miracle” has been interpreted much later as stemming from the 
relabeling invariance of the variational formulation of the Euler equation and use of Noether’s theorem (Frisch & Villone 2014, Section 5.2). 
The resulting Cauchy invariants equations have normally the initial vorticity in its r.h.s but, in the present case, the vorticity vanishes and the 
equations take the following form: 

(i = l,2,3) (14) 


where Sijk is the fundamental antisymmetric tensor, and summation over repeated indices is assumed. The l.h.s. of d is a Lagrangian curl. 
The second part makes use of all three equations ( I 12 ab -( ll 2 cb . and gives a third scalar equation: 


SiklSjmn^ nX I 


(l + Aa®) dja + 




(15) 


This is done by taking the Eulerian divergence of l ll2ab to obtain the Eulerian Laplacian of the gravitational potential, which is then substituted 
in l fT2cb . The density contrast <5 is expressed in terms of the Jacobian, using l ll2bb and also substituted in the r.h.s. of il2cb . Finally, the 
Eulerian divergence is expressed in terms of Lagrangian space derivatives and of the inverse Jacobian matrix. The latter is given in terms 
of the Jacobian matrix using the following identity: = EikiEjm-n^^xP^^xi/{2J), which follows from the observation that the 

Jacobian matrix of the inverse Lagrangian map is the matrix inverse of the Jacobian of the direct map (see, e.g., Buchert & Goetz 1987). This 
derivation fails if the Jacobian J is not everywhere strictly positive in the space and time domain under consideration. In the latter case, the 
Lagrangian formulation d-d becomes invalid. We shall come back to this issue in Section lh^ 


4 THE FORMAL TAYLOR EXPANSION AND THE RECURSION RELATIONS 


In the LPT literature, analytical derivations of results for ACDM usually rely on an expansion in powers of small displacements, which 
amounts to performing a power-series expansion in the initial peculiar gravitational potential (see e.g., Matsubara 1995; Bemardeau et 
al. 2002). In this paper we follow a different approach and perform an expansion in powers of a-time. Note that for the case of an EdS 
universe, one can obtain recursion relations for the displacement by both expansion techniques, and the resulting recursion relations are 
formally identical (Rampf 2012). For a ACDM Universe, however, explicit results for the displacement will generally differ depending on 
what expansion scheme is used. 

We thus seek a solution to the Lagrangian equations d-d in the form of a power series in the o-time (i.e., a Taylor series) for the 
displacement ^ = x — q, 

OO 

i{Q,a) = (16) 

S=1 

At first order, we have 

= (17) 


as follows from v = d^x and d- The convergence of the series d will be examined in the next section. The Jacobian of the Lagrangian 
map J = det(/ + where I is the identity matrix, can be written as a sum of four terms. After substituting the expansion H 6 b . the 

Jacobian becomes 


J=l + ^^Wa'’+ ^ ^ 

s = l ni-\-n 2 ~s ni+n2+7i3 = s 


(18) 


where the sums are restricted to values of ni, n 2 and 713 that are strictly positive and where the various quantities /ri, /r 2 and /is are 
expressions linear, quadratic and cubic, respectively, in the Lagrangian gradients of the Taylor coefficients, given by 


fiY = 

(19) 


( 20 ) 

,Xni,n2,n^) _ 1. . yyL ^(ni) YyL^(n2) Y7L^(n3) 

^3 — Q^tkl^jmn y m^k ^ n^l ^ j 

( 21 ) 


Substituting the expansion ( 116b into ( 115b and identifying the coefficients of the various powers of a, we find, for s > 0, 


s + ^ ) (s - 1) - 


{ [f - 2n2 (n 2 + 0] - 2A(n2 - 3) (n 2 - 1)| 


- A{s - 3) {s - 1) fj.Y 


ni+n2=s 

(s-3) 


ni+n2+n3=s 


- — 3n3 ns + - 




( 22 ) 


Here, by construction, all the /i coefficients vanish if one or several of their upper indices are zero or negative. We symmetrise the r.h.s. and 
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obtain a sequence of relations for s > 1: 

>("S) _ (•s) vtL (init) cs \ ^ ^ — 3) 


+ E 

0<n<s 


(3 - s)/2 - - (s - n)^ ,,(n,s-n) , ^ 4s - - (s - n)^ - 6 (n,s-n-3) 

M 2 + A , n lc)\( 1 \ M 2 


(s + 3/2)(s- 1) 

2 _2 _2 


(s + 3/2)(s-l) 

^2 I _2 I _2 


(3 s)/2 72]^ 722 723 (tt-i ,Ti2 1^3) A ^2 ~l~ ^3 ^ ('^11‘^2 ~3) 

(s + 3/2)(s- 1) “ (s + 3/2)(s-l) 


+ E 

'n.i+n,2+n3 = s 

where Sf is the Kronecker symbol. Similarly, we obtain from the Cauchy invariants equation ifTit . for s > 1 


V X = T 


(s) — rri(s) _ 1 


2 ^ >, 
l<fe<3, 

0<n<s 


^ E 


( 23 ) 


(24) 


As we see, (123} and (IS} give the Lagrangian divergence and curl of the a-time-Taylor coefficient of order s in terms of lesser-order ones. 
We observe that the first explicit occurence of A in the recursion relations is at the fourth order, which is the same order at which the linear 
growth D{a), at short a-times, differs from its EdS value: D(a) = a — (2/ll)Aa'* + 0{U). 

Eqs. ( 123b and (124b give a Helmholtz-Hodge decomposition of into curl and divergence, from which it can be retrieved in a standard 
way (cf. e.g. Arfken & Weber 2005), namely 


^(S) ^ ^-2 _ ,^L ^ 


(25) 


where is the inverse Laplacian in Lagrangian coordinates (with the boundary conditions taken into account), and and denote 
the r.h.s.’s of Eq. ( 123b and Eq. ( 124b . respectively. 

Using ( 125b and considering the first-order Lagrangian derivatives of we can combine the two recursion relations ( 123b and (|24} into 
a single recursion relation for the gradient tensors of the time-Taylor coefficients, which reads 




s-3 
s-f 3/2 


e E 


l<j<3. 


( 



<1 

1 

\ l<fe<3, 

\0<n<s 



( 


+ Cii 


E 

0<n<s 


(3 ^)/2 72 (s 72) (n,s —ri) , a ^ ^ (n,s —ri — 3) 

(s-b3/2)(s-l) ^ (s-f3/2)(s-l) 


+ E 

^l+^2+^3 = ^ 


/\/222 222 

(3 Sj/2 72^ 722 ^3 (^1)'^27^3) I ^ ^2 ^3 ('n.i ,n 2,713 —3) 

(s-b3/2)(s-l) ^ (s + 3/2)(s-l) 


(26) 


/ 


where Cij = is an operator of the Calderon-Zygmund type (cf. Zheligovsky & Erisch 2014). The recursion relations d26b are 

explicit but look somewhat lengthy. Actually, they enjoy an important property which allows us to prove time-analyticity in the next section: 
for s > 1, it is easily checked that all the coefficients involving s, n\, ... are bounded (rational) functions. In those not involving A, the 
coefficients are bounded by unity and the other ones by A > 0. 


5 THE TIME-ANALYTICITY OF THE LAGRANGIAN SOLUTION 

After having found explicit recursion relations for the time-Taylor series of the displacement field, we are naturally led to ask: is 

CX3 

= (27) 

S=1 

a convergent series that thus defines a time-analytic function in the neighbourhood of the origin? 

Let us give first some overall ideas on how the convergence results will be established. Given that the recursion relations take their most 
compact form ( 126b in terms of gradients, it is simpler to first prove the convergence for the gradient series of 

OO 

= (28) 
S=1 

From the recursion relations ( 126b for the gradients of the Taylor coefficients, we shall be able to derive polynomial recursion inequalities 
for their norms. By reintroducing the o-time and summing over all orders (i.e., by using a generating function), we shall obtain a single 
inequality for an o-dependent cubic polynomial. The study of the evolution of its roots will give us the required bounds on the norms of the 
time-Taylor coefficients of the gradient series, from which the convergence of the series ( 128b and thus time-analyticity are established. The 
analyticity of the Lagrangian map and of the displacement field are then a consequence. 
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Let us proceed now with the details. First, we observe that the r.h.s.’s of the recursion relations l l26t for the gradient tensors are them¬ 
selves mostly polynomials (of degree not exceeding three) of the lesser-order gradient tensors. We write “mostly” because there are also 
the Calderon-Zygmund operators, Cij. In the present case, these operators stem from the nonlocal nature of the gravitational interaction. 
Technically speaking, these are pseudo-differential operators of degree zero (because an inverse Laplacian is compensated by two space 
derivatives). Essentially, such operators do not change the degree of differentiability of the functions they are applied to. But this is true only 
when using a suitable function space in which the Calderon-Zygmund operators are bounded. 

For our purposes the function space in which such matters are simplest is the space of spatially periodic functions, say, of periodicity 
27r in all three space variables and such that their Fourier series is absolutely summable. Concretely, let a periodic function f(q) (of scalar, 
vector or tensor type) be expanded in a Fourier series: 

f(q) (29) 

k 

where the summation is over all triplets k of signed integers. A periodic function is said to be in ^i, if the norm 

l|f||=I]|fte| (30) 

k 

is finite (for vector and tensor quantities, the modulus is defined as the maximum over all indices of the modulus of the various components). 
Since, for any i and j, one has \kikj/k‘^\ < 1, the Calderon-Zygmund operators are obviously bounded by unity in this space. It is also 
elementary to show that this space enjoys the algebra property: 

l|fg|| < l|f||||g||, (31) 

for any pair of functions f and g of finite norm. Other function spaces with similar properties may be more realistic for cosmological 
applications in so far as periodicity is not required, such as Holder spaces. For this we refer the reader to Section 2.5 of Zheligovsky & Frisch 
(2014). 

Given the structure of the recursion relations and using the boundedness of the Calderon-Zygmund operators and the algebra 
properties of the £i norm, it is elementary to show that if is in the space £i, so will be the gradients of the Taylor coefficients 

jjjg assumption that the gradient of the initial velocity has absolutely summable Fourier coefficients is the key hypothesis of the 
present work. 

Knowing that all the (Lagrangian) gradients of the Taylor coefficients are 'm£i is not enough to ensure the convergence of the gradient 
Taylor series. For this, we need to obtain suitable bounds on the £\ norms of and then we need to show that these bounds imply 

the convergence of the gradient time-Taylor series. Here we work with the time-Taylor series for the spatial gradient of the displacement and 
establish its time-analyticity, from which follows the time-analyticity of the displacement (and also, of course, of the full Lagrangian map). 
To obtain the bounds on the £i norms, we shall mostly follow the approach presented in Zheligovsky & Frisch (2014). Of course, due to the 
presence of the cosmological constant, the third-order polynomials contain additional terms, and their study is more involved. 

The first step is to bound the various || ||, using the boundedness of the Calderon-Zygmund operators, the algebra property, and 

the boundedness of the rational coefficients, as explained at the end of section|4] We use d26b and thereby obtain, for any s > 1: 

-l-3A||v^^^““®^j| 4-12 ^ ||j| ^ ^ 

+ ® + 6A ^ . (32) 

i+j-\-k = s i-\-j-\-k = s 

The second step is to introduce a generating function using gradients of Taylor coefficients, namely 

00 

C(a) = ^||v^^("^||a“. (33) 

S = 1 

It is easily checked that this generating function is an upper bound for the £\ (absolutely summable Fourier series) norm of the gradient of 
the displacement field at time a. Multiplying l l32b by 0 “ and summing over s from one to infinity, we obtain the following inequality for the 
generating function: 

c < a -h 12C^ -f 6C® + 6Aa®C^ + 6Aa®C^ -f 3Aa®C , (34) 

where ^ stands for C(a)- To study this cubic polynomial inequality, it is convenient to introduce rescaled variables defined as follows: 
as , A s , C(o) = C(a) • (35) 

The polynomial inequality d34b becomes then 

PA (a, C) = 6 (1 -f \c?) C® + 6 (2 -h An®) C" + (3Ao® - l) C + a > 0 . (36) 

Sketches of the cubic polynomial PA(a, for A = 0 and for A yi: 0 are respectively given in Figs .Q] and |2] 

The cubic polynomial pa (a, has at least one real root. For o = 0 it has three roots: zero, a positive root and a negative root. Since 
the displacement, and thus C,, have to vanish at time a = 0, the only physically relevant root of the cubic polynomial is zero. Now, suppose 
we let the (rescaled scale-factor) time 0 be small and positive, then the physically relevant root moves from zero to a small positive value 
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Figure 1. Sketch of the polynomial (36), for A = 0, i.e., for the case of an EdS universe. Shown are three values of the rescaled time a for which 0 < 
a< < T(a = o) < a> , illustrating the behaviour of real roots of po(cij C) (the roots are shown as the points of intersection of the graph of po(ci< 5 C) ^ttd the 
horizontal axis). On increasing a, the graph slides up as a rigid curve. As a result, the roots (^i and ^3 move to the left (i.e. become smaller), whereas (^2 moves 
to the right (i.e. becomes larger). At the critical value T(0), when dpo/dC, = 0 and the discriminant A vanishes, the two roots <^2 ^J^d <^3 collide and then 
disappear (with the emergence of a pair of complex conjugate roots). 


C 2 (ci) ~ a. The polynomial p\{C) is then strictly positive over the whole open interval ] 0 , C 2 (ti)[, as required from the inequality < [3^ . and 
this is the only physically relevant branch of positivity. Hence, within the interval [0, a], the values of the generating function ^ are bounded 
by C 2 (ti). As we increase the value of a, this boundedness property will hold as long as there are two positive roots in for the polynomial 
Px{a, C). This is true until the two positive roots merge into a double root and then turn into complex roots. The value of the rescaled time 
for which l l36t has a double root are the zeros of its discriminant, i.e., 

A(o, A) = 12 (14 - 684n - 81o^ - 76a®A - 702n'‘A - 162o'^A + 75o®A^ - 81o®A^ + 90o®A® + 90 o“A® - 27o^^A'‘) = 0 . (37) 

Because A(0, A) > 0 and the highest-order term, proportional to has a negative coefficient, it immediately follows by continuity that the 
discriminant equation has at least one positive root. For our purposes we need the smallest of such positive roots, since it is an upper bound 
for the physically relevant branch of pa(C) > 0 : 

r(A) = inf Ui I A(ai,A) = 0, (38) 

i 

called here the critical time. (In Section lfidl we shall show how T{\) can be calculated.) 

The importance of the critical time lies in the following result:/or any |n| < T(A) the time-Taylor series for the gradient of the 
displacement (now expressed in terms of the rescaled time) is convergent. (Here |a| denotes the modulus of the rescaled time a which can 
now take also complex values.) Indeed, it follows from <331 and <351 that 

OO 

= (39) 

S = 1 

where = || ||/|| ||® > 0. IfC(r) is bounded by some constant C, we have 

OO 

C(r) = y^c;(®)T® < c. 

s=l 


( 40 ) 
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Figure 2. Sketch of the polynomial <36t for any fixed rescaled time variable a > 0 and A 7^ 0 (i.e., the ACDM case). Three values of lambda are shown for 
which A< < A* < A>. The value A* denotes the critical value for which the cubic polynomial px develops a double root, i.e., where (^ 2 (^(-^)) = C3(^(^))- 
In that case the fixed value of the rescaled time a is precisely the critial value T(A*) for which the generating function of the displacement is bounded. On 
increasing A for any fixed rescaled time a > 0, the graph gets shifted into the upper left direction. 


Since all the terms in the sum are non-negative, it follows that 

^(s) ^ ^ 

and thus 

00 

(42) 

S = 1 

For |a| < T the Taylor series HU is hounded by a convergent geometric series. We have thus shown the convergence of the Taylor series in 
the complex time plane inside a disk, centered at the origin, o = 0 of radius at least T{\). The actual radius of convergence of the Taylor 
series, called the radius of analyticity, for which T(A) is only a lower bound, in general, can only be determined numerically. 


6 FURTHER RESULTS ON ANALYTICTY 

6.1 The lower bound on the radius of analyticity and its dependence on the cosmological constant 

The lowest positive root T{X) of the discriminant equation ( 1371 gives us a lower bound on the radius of analyticity of the time-Taylor 
series for the Lagrangian map. This discriminant equation is of twelfth degree in the rescaled time variable a = a|| and we have 

not been able to solve it explicitly by radicals. We can however solve the discriminant equation perturbatively for small and large rescaled 
cosmological constant A = A/|| ||® and numerically for other values. 

For A = 0, we recover the results for an EdS universe (Zheligovsky & Frisch 2014). In that case the discriminant equation reduces to 
A(o, 0) = 12(14 — 684a — 81n^) = 0, whose lowest positive root is 

yEds _ ^ 3^2 _ 38/9 ~ 0.0204 . 


(43) 
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T{\) 



A 

Figure 3. Numerical solution of the discriminant equation . The results agree with the small- and large-A expansion EH and EiJ in the appropriate ranges. 


For small A, it is then straightforward to obtain the A-dependence of this root 

T(A) = 3^~_^^g^yEds ^ ^ 0.0204 - 1.1197 • 10“® A + 0(A") . (44) 

For large positive A, we make the changes of variable a = aA^^® and A(a, A) = A(o, A)/12, and rewrite the discriminant equation as 
A(a,A) = P4(a®) + A-^/®Pio(a) + A-"/®P8(a) =0, (45) 

where 

P 4 (a®) = 14-76n® + 75a®+90a’^-27o^^ Pio(a) =-684a - 702a'‘+ 90a“ , P8(a) =-815^ - 162o'^ - 81a® . (46) 

The polynomial P^icd) has a positive double root at Oc = from which follows that, to dominant order, r(A) is proportional to 

The first subdominant correction is obtained by Taylor expanding the polynomial P4. Pro and P8 around he. This eventually 
reveals the following large-A behaviour: 

r(A) = 3“^/®A“^/® + 0(A“^/^). (47) 

We have solved the discriminant equation ( 137b numerically for more than one hundred values of A, suitably distributed between 10“^ 
and 10^^. The results are shown in Fig.[^and agree with the small- and large-A expansion ll44b and d47b in the appropriate ranges. Note that 
r(A) is a monotonically decreasing function of A. 

Finally, we observe that, within the functional framework of the space £i of absolutely summable Fourier series, the precise values of 
the constants appearing in the low-A expansion and the large-A expansion (113 can most likely be improved. This will result in longer 
time intervals of guaranteed analyticity. For this one should adapt to the cosmological context the detailed Fourier-space estimates found, for 
incompressible flow, in Section 2.3 of Zheligovsky & Frisch (2014). 
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6.2 Shell-crossing and analyticity 


In section [3 we pointed out that the Lagrangian formulation and the subsequent time-Taylor expansions as carried out in the present paper, 
are invalid if the Jacobian J vanishes during the interval of guaranteed analyticity. The vanishing of the Jacobian corresponds to the crossing 
of various fluid-particle trajectories. It is generally known as “shell-crossing”. We now show that during the time interval of guaranteed 
analyticity, no shell-crossing can occur. Actually, denoting as usual the Jacobian of the Lagrangian map by J, we shall show that | J — 1| < 
—37/9 -I- 3x72 « 0.132, which trivially implies the non-vanishing of the Jacobian (and not even a close call). 

For this, we observe that the Jacobian can be written as follows in terms of the displacement field 

J = det(/ + V'^0 = l + V'^-^-F +det(V^O, (48) 

where I denotes the identity matrix. Using the Taylor expansion of the displacement in powers of the a-time and the generating function 
/ = I we obtain the following bound 

|J-1| <6C®+6C"-b3C = 0(0- (49) 


For convenience we have plotted the cubic polynomial Q{Q in Fig.|4l 

As we have seen in Section|5] for any given rescaled cosmological constant A and rescaled time n, the permitted values of the generating 
function / are between zero and the smallest positive root 1 / 2 ( 0 ) of the cubic polynomial Pa(C, q), given by i36t . We claim that this root is 
an increasing function of a for fixed A. Indeed, Pa(q, ( 2 ( 0 )) = 0, over the whole range of relevant o values. By differentiating it with respect 
to 0 , we obtain: 


9 pA(n,C 2 (a)) _ 9/2(0) apA(a,C) I 


-t- 1 -f 3Ao ( 6/2 T 6/2 “b 3 / 2 ) — 0 . 


9o 9o d( lC=C 2 (c) 

Since [1 -|- 3Ao^(6/| + 6/| + 3 / 2 )] > 0 for (2 > 0, we infer that 
9/2(0) 9 pa(o,/) I 


9o 


9 / 


k=C2(a) 


<0. 


(50) 


(51) 


The positivity of follows from < 0. The latter is a consequence of the fact that the polynomial pa(o, /) has three 

“ k=C2(a) _ 

real roots, for fixed A and fixed o in the range of analyticity. Inspection of yb) shows that their product is negative, that one (/i) is negative 
and that two (/2 and (3 ) are positive. The derivative is obviously positive for large / > 0 and changes sign somewhere between /2 

and (3. Since (2 < /a, this derivative is negative at / 2 . This proves the monotonic increase of (2 with q. 

Next, we use to bound IJ — 1| by 6/1 + 6/1 + 8 / 2 . The largest possible value of (2 is obtained at the critical value a = T(A) 
where the polynomial px(a, /) has a double root in /. This value can be obtained by demanding the simultaneous vanishing of pa(ci, /) and 
of its derivative with respect to /, i.e., Pa(/ 2 ) = Pa(/ 2 ) = 0 and (2 > 0, where a prime denotes a partial derivative with respect to /. From 
these conditions we obtain 


^ C,(T(A)) ^ . (52, 

Since this is a decreasing function of A, it is bounded by its EdS value (A = 0), where the latter is /2 = — | -b |\/T 8 ~ 0.04044. It follows 
that IJ — 1| < 0.132, which is well below the value unity where shell-crossing might occur (see the vertical thick line in Fig.|4j. 


6.3 Analyticity in the linear growth function D 

As already stated in Section[2l it could be of interest in cosmological studies to use as time variable, not the cosmic scale factor a, but the 
linear growth function D(a), which is the growing solution of ([7]l. With the normalisation lima->o D{a)/a = 1, the linear growth function 
is given by 

D{a) = ax/l + Aa^2F3 Q, y, -Aa') . (53) 

where 2 F 1 is the Gauss hypergeometric function. This solution is analytic around a = 0 and has the small-a expansion D{a) = a — 
(2/11) Ao^ -b 0(U). (The full expression can be obtained from the Taylor expansion of the Gauss hypergeometric function, cf. Gradshteyn 
& Ryzhik (1965), pp. 1039-1045.) As D{a) is analytic around a = 0 and also its derivative dD jda does not vanish at a = 0, it follows 
that D{a) is invertible at around a = 0. The inverse linear growth function a(D) is also analytic and has the following low-order expansion: 
a(p) = D + (2/11) + Oip^). 

We have shown in Section|5]that the Lagrangian map a;(q, a) is an analytic function of a, near a = 0. The composition a:(q, a{D)) of 
the Lagrangian map with the inverse linear growth function is also analytic in D. (This follows basically from the observation that analytic 
functions are complex-differentiable functions and the use of the chain rule.) We remark that this argument is not valid when the cosmic time 
is chosen as time variable, instead of D, as the relation between a and t is not analytic near the origin. 

Can we obtain simple recursion relations for the ZJ-Taylor coefficients of the displacement /, as we have done in Section|4]? Recursion 
relations, yes, but simple, no. Indeed, from d-GIll, we can derive a set of Lagrangian Euler-Poisson equations in the time variable D. 
However the D-dependence is then essentially not polynomial, but will involve the full inverse linear growth function a{D), whose Taylor 
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0(C) 



c 

Figure 4. Lower bound on the norm of | J — 1| < Q(C) = + 3C, see Eq. The vertical thick line at roughly ^ = 0.04044 indicates the largest 

possible value of the smallest positive root C, 2 {T{x = o)) of the polynomial C) which analyticity is guaranteed. Since (^ 2 (T'(A)) is a decreasing 

function of A (see (52)), it follows that | J — 1| < 0.132 for any value of A > 0, which is well below the value unity where shell-crossing might occur (the 
dotted line). 


expansion has an infinite number of terms. As a consequence, the recursion relations will be very involved. If one really wants to reexpress 
the displacement as a power series in D, it is much simpler to just substitute the Taylor series for a{D), obtained by inversion of the Taylor 
series of D{a) into the expansion jl6b in powers of a. Note that, with the latter strategy, there is no need to solve a succession of time 
differential equations for the nth order growth function (usually denoted with E, F ,.because the actual expansion parameter is precisely 
the cosmic scale factor, used as a time variable. 


6.4 Effect on analyticity of the inclusion in the Friedmann equation of curvature and radiation terms 

Inflation theory and observational evidence since the late nineties favour a flat Universe with zero curvature [see, e.g., Linde 1984 and 
Netterfield et al. (Boomerang Collaboration) 2002]. It is nevertheless of interest to point out that our analyticity result for a ACDM Universe 
survives when a small amount of curvature is introduced in the Friedmann equation (H. The latter then becomes 

=a“^ + fca“^ + A, (54) 

where k is here not limited to being 0, ±1 (although according to Mukhanov (2005), the only trivial background geometries which do not 
spoil exact homogeneity and isotropy are indeed the ones with curvature 0, ±1). We have repeated the analysis of Sections |2)[5] with this 
modified Friedmann equation and found that very little is changed. The slaving conditions ( IlOb are unchanged. The recursion relations ( 126b 
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are minimally modified. The polynomial inequality ( 136b in the generating function C, takes now the following form in rescaled variables: 

c® + 6 ^2 + An® + 

where R — |fc|/|| and o = a|| and A = A/|| as before. Analyticity in the (rescaled) cosmic scale factor 

a is proved as before, by studying the discriminant equation associated to the cubic polynomial C)- Indeed, this is again an equation 

of twelfth degree in a, whose lowest and highest degree terms, unchanged from the ACDM case, are respectively positive and negative, thus 
implying the existence of at least one positive root. 

All this is not surprising; Actually, the curvature term in the Friedmann equation involves the inverse second power of the cosmic scale 
factor, which for a —>■ 0 is subdominant with respect to the matter term and can thus not bring about much change. Very different 
would be the inclusion in the background of radiation ejfects, stemming from a black-body distribution of photons. This adds a term in 
the Friedmann equation, proportional to the inverse fourth power of the scale factor, which drastically changes everything: slaving is lost, 
analyticity around a = 0 is lost, etc. The current favoured ACDM model has a background radiation density well below the matter density 
at decoupling and can actually be discarded after decoupling. When relating our results to numerical results or observations, it may be more 
convenient to use the cosmic time t rather than the a-time. Of course, when performing that change of variable, the radiation term in the 
Friedmann equation does not need to be omitted. 


C) = 6 (l + + 2 .ffa) 


+ (sAo^ + yjfa-1 ) C +0 > 0, 


(55) 


7 CONCLUDING REMARK S 

In a Eulerian framework space and time are intertwined by Galilean transformations (in the Newtonian approximation). It follows that, with 
the limited spatial smoothness assumed in the present work - roughly, that the initial density fluctuations are slighly better than continuous 
in the spatial variables - one cannot obtain better than a corresponding limited temporal smoothness. In a Lagrangian framework, as we have 
seen, the situation is drastically different: provided we use an appropriate time variable, such as the cosmic scale factor a or the linear growth 
function D. Obtaining analyticity in time of the Lagrangian trajectories requires then only a limited initial spatial smoothness. 

Even more surprising is that this analyticity is around the point a = 0, which seems to extend this utterly smooth birth of our structured 
Universe to the very origin of time, well before decoupling of matter from photons. Of course, this is just a property of the mathematical 
model used here to describe structure formation, namely the Euler-Poisson equations in a ACDM Universe. The interesting feature is that this 
model allows us to have a well-posed problem, devoid of catastrophic behaviour near decoupling, provided we use the slaving conditions (Uni 
at a = 0. These slaving conditions resemble those used to start W-body simulations at early times, (i) small density fluctuations, (ii) an initial 
curl-free peculiar velocity field, related to the gravitational potential as in the ZeTdovich approximation. That these slaving conditions can 
be extended all the way to a = 0 means, from the point of view of boundary layer analysis (Cole 1968), that the matter-dominated era 
constitutes an outer solution for which, to leading order, the boundary matching to the inner solution (describing the primordial era) can be 
replaced by just an initial condition. This is true as long as, in the Eriedmann equation for the background evolution in the matter-dominated 
era, we include only a matter term, a cosmological term and a possible curvature term, but no radiation term. 

We remind the reader that our way of proving analyticity in the a-time is fully constructive. It rests indeed on a set of novel all-order 
recursion relations l l26t for the time-Taylor coefficients of the displacement of (fluid) particles. Not only are these relations fully explicit, but 
they have a very specific mathematical structure, allowing us to obtain bounds for all the Taylor-coefficients and thus to establish analyticity 
by elementary means. 

We now briefly indicate some of the possible future developments exploiting the analyticity result and the recursion relations. We shall 
mention two interlinked topics. One is about numerical integration of the Euler-Poisson equations by a multi-step technique and the other 
one about cosmological reconstruction. 

The Taylor method as described in Sections |4] and [fallows us to determine the solution of the Euler-Poisson equations for any time 
0 < a < i?(0), where R{0) is the radius of convergence of the Taylor series around a = 0. (Here, it is best not to work with the rescaled 
time a.) A procedure, inspired from the Weierstrass analytic continuation technique for analytic functions, may allow us to obtain the solution 
beyond the time i?(0). Eor this, we select a sequence of times 0 < ai < 02 < as < ... in such a way that a„+i is within the disk of 
convergence of the time-Taylor series around the time a„, that is such that |a„+i — a„| < R{an), where R{a) is the radius of convergence 
of the Taylor expansion of the solution around time a. What we here call the solution is not just the Lagrangian map, but includes also the 
density which is controlled by the inverse of the Jacobian. Indeed, if the Jacobian vanishes, shell-crossing takes place and the Lagrangian 
map ceases to be invertible. Invertibility is essential, because at the end of each step, we must revert to a Eulerian description to be able to 
make a fresh start. Methods combining a Lagrangian step with a reversion to Eulerian coordinates are called semi-Lagrangian. A detailed 
description of a semi-Lagrangian method, called the Cauchy-Lagrangian method, which makes use of Cauchy’s Lagrangian formulation of 
the ideal fluid dynamics can be found in Podvigina et al. (2015). 

When dealing with the cosmological Euler-Poisson equations, the method of Podvigina et al. (2015) must be suitably adapted. One of 
the new difficulties stem from the non-autonomous character of our Euler-Poisson equations: the time appears explicitly in the equation. As 
a consequence, the Euler-Poisson equations are not invariant under a-time translations, but this problem can be circumvented. 

In so far as the regime described by the Euler-Poisson equations does not include multi-streaming, a Cauchy-Lagrangian numerical 
method will not be able to address the same questions as W-body simulation techniques. There is however one important problem for which 
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it appears well suited, namely the reconstruction of the dynamical history of the Universe from present-epoch galactic surveys (Frisch et 
al. 2002). It has heen shown hy Brenier et al. (2003) that the solution to the Euler-Poisson equations is uniquely determined hy two boundary 
conditions, the slaving conditions 01 at the initial time and the density of matter at the current epoch. The latter is obtained from large 
galactic surveys. Euler-Poisson reconstruction is an extension of the variational A'^-body reconstruction, introduced by Peebles (1989) for 
handling galaxy data on relatively small scales, such as that of the Local Group. Euler-Poisson reconstruction, which is meant for significantly 
larger scales where a continuum description is applicable and where multi-streaming can be ignored, does not suffer from the non-uniqueness 
problem of the variational A^-body reconstruction. So far, Euler-Poisson reconstruction has not used the full solutions of the Euler-Poisson 
equations but just low-order Lagrangian perturbation approximations, such as the Zel’dovich approximation or the next order, denoted by 
2LPT. Within such approximate frameworks, reconstruction becomes a Monge optimal transport problem (Erisch et al. 2002) with quadratic 
cost, which, after discretisation, can be solved by very efficient assignment algorithms such as the auction method of Bertsekas (1992). 
Alternatively, one can use Brenier’s theorem (Brenier 1987) and solve an equivalent 3-D Monge-Ampere equation non-linear PDE by 
iterative techniques (Zheligovsky, Podvigina & Frisch 2010). Reconstruction by such methods is known as Monge-Ampere-Kantorovich 
(MAK). We propose to use the MAK reconstruction as a starting point of full Euler-Poisson reconstruction, applying an iterative Newton- 
type method in which at each stage a Euler-Poisson initial value problem is solved by a Cauchy-Lagrangian method. 

Einally, we remind the reader that our proofs were given entirely within a Newtonian framework. In that framework the time-analyticity 
results revealed a strong asymmetry between space and time in Lagrangian coordinates. The reason for that strong asymmetry is the “missing” 
convective term in the Newtonian fluid equations, when formulated in Lagrangian coordinates. (In the notation of our Eulerian approach in 
Section|2 the convective term is v-'V^v.) An interesting question is whether our time-analyticity results could be generalised to a Lagrangian 
coordinates formulation of a curl-free dust-fluid in General Relativity. It is well-known that the so-called synchronous-comoving coordinate 
system corresponds to a relativistic Lagrangian frame of reference, and, similar to the Lagrangian coordinates approach in the Newtonian 
framework, there is indeed no convective term in the relativistic equations of motion (in the so-called ADM approach; see e.g. Matarrese & 
Terranova 1996). These relativistic fluid equations are however decorated with Ricci-tensors that involve (single and double) spatial gradients 
acting on the metric (and on its inverse), and this could well imply a limited temporal smoothness for limited spatial smoothness. Besides that 
unsolved problem, to obtain a time-analyticity result in General Relativity along lines similar to those of the present paper, it seems necessary 
to first obtain explicit all-order recursion relations. A possible starting point for such investigations could be the Lagrangian equations recently 
obtained by Alles et al. (2015). 
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APPENDIX A: DERIVATION OE THE PECULIAR EULER-POISSON EQUATIONS WITH COSMOLOGICAL CONSTANT 


Here, although entirely elementary, we highlight the derivation of the peculiar Euler-Poisson equations for the sake of completeness. 

We begin with the (full) Euler-Poisson equations i fTalfT^ . which we repeat here for convenience: 


dtU + 


dtQ + V 


iy-Vr) 
r ■ (qU) = 


U = -Vr 


= dyrGp — 3A . 


(Al) 

(A2) 

(A3) 


Here U = U{r, t) denotes the velocity field, (j)g = (j>g{r, t) denotes the cosmological potential, q = g{r, t) is the fluid density, and here the 
cosmological constant is 3A. Let us solve these equations for the background variables in an exactly homogeneous and isotropic Universe, 
where q = g{t), U = H{t)r = (a/a)r, and (pg = (j)g. We substitute these expressions into the above equations. For the divergence of 
EQ.dATTl. we then have 

sA + = -Vl4>g , (A4) 

whereas for Eq. dA3b . we have at the background level 

Vl^g = TttG^ - 3A . (A5) 

By substituting the last expression into Eq. ( IA4b . we obtain 

a AtvG . . ,,,, 

- =- 5 —P +A- (A 6 ) 

a 3 

This is the usual form of the second Friedmann equation for a pressureless fluid. In the units after Eq.l[3i in the main text, the above is 
a/a = —a~^/2 + A, whose solution is a{t) = A“^^®[sinh(3'\/Af/2)]^^®. 

Now we take the fluctuations of our variables into account, i.e., we demand 
a{t) 


g = g{t) [1 + 5 ] , 


a{t) 


(A7) 


9g ^ 9g + 'fig ^ 

In order to write the Poisson equation for the peculiar potential pg, we use Eqs. ( IA3b and ( IA5b . and obtain 

'V^ifg = = TttG (J) — p) = . (A8) 

Observe that A has dropped out in the peculiar Poisson equation, i.e., A does only affect the background evolution of the Poisson equation. 
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Let us now derive the peculiar evolution equation of Eq. jAU by use of JA7b . We have 


a. 


(% + an) 


+ 


-r + aii I • Vt 
a 


— r + au I = — Vr 
a 


+ ‘fig] ■ 


(A9) 


Let us do the calculations in JA9I > step by step. The l.h.s. of JA9I > can be written < 


l.h.s. = —r- -r + au + adtu H— ttU ■ Vt-It- + a(u ■ Vrlf + a(r ■ ^r)u + (u ■ ^r)u 

a a-* ' ' ' ' ' ' ' ' 


... 2 -2 

Uj , -w ^ -w 61 * ^ \ 2 / \ 

= —r- -r + au + aotu H —-r + au + a(r ■ Vrlu + a (u ■ Vrju . 

a a-‘ a-‘ \ \ / 

We thus obtain 


(AlO) 

(All) 


adtu + a [u ■ \^rju + 2au + d(r • Vr-jit = —'Vr'fig- 
Now, we change the dependence of the above functions from r to a; = r/a. For an arbitrary function f{t, r) = f{t, a{t)x), we have 
dtf\xfixed ~ dtf\r fixed a(^X ■ f\tfixed • (A12) 

Doing the same for our Eqs. dAl lb . and note that Vt- = ^ Vx, we obtain 

1 1 

adtulrBxed — air ■ Vr)u + a^iu ■ -VFju + 2du + dir ■ Vr)u =- Vx<Po- (A13) 

a a 

Then, we divide by a and obtain 

dtU + lu ■ 'Va:)u + 2-U = --^Veefig ■ (A14) 

The above equation, for an EdS universe was derived by Brenier et al. (2003). Note that we have corrected a typo in their equation (A12). 






